home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / HyperCard Files / MacMathPascal / card_3829.txt < prev    next >
Encoding:
Text File  |  1988-09-05  |  6.9 KB  |  239 lines
  1. -- card: 3829 from stack: in
  2. -- bmap block id: 0
  3. -- flags: 0000
  4. -- background id: 2607
  5. -- name: ALGAG3D.PAS
  6.  
  7.  
  8. -- part contents for background part 22
  9. ----- text -----
  10. ALGAG3D.PAS
  11.  
  12. -- part contents for background part 23
  13. ----- text -----
  14. ALGAG3D.PAS
  15.  
  16. -- part contents for background part 17
  17. ----- text -----
  18. ALAG3D, PLOT3[function]
  19.  
  20.  
  21.  
  22. PURPOSE 
  23.  
  24. ALAG3D performs an interpolation of the function U=U(x,y) defined on a 4 x 4  or larger grid by performing a weighted averaging of eight Lagrangian interpolations in the area of the point of interest. The X values and the Y values need not be  uniformly spaced and ΓêåX need not equal ΓêåY.   ALAG3D  requires the passage of an X and Y array, the  X and Y array size, a 2d Z array and the points XIN, YIN for which the interpolation is to be done.
  25.  
  26.  
  27. TYPE REQUIREMENTS
  28.  
  29. CONST
  30. NUMX = maximum number of X data points-1;
  31. {may be larger than XSIZE}
  32. NUMY = maximum number of Y data points-1;
  33. {may be larger than YSIZE}
  34.  
  35. Type
  36. FLOAT = REAL or DOUBLE or EXTENDED;
  37. DATAX = ARRAY [1..NUMX] of FLOAT;
  38. DATAY = ARRAY [1..NUMY] of FLOAT; 
  39. RECTANG = ARRAY [0..NUMX,0..NUMY] of FLOAT;
  40. FIXXED = INTEGER or LONGINT;
  41.  
  42.  
  43. CALLING PROCEDURE
  44. var
  45. X:DATAX;{Known X values}
  46. Y:DATAY;{Known Y values}
  47. Z:RECTANG;{Known Z values}
  48. XSIZE,YSIZE FIXXED;{Size of X AND Y array}
  49. X0,Y0,Z0:FLOAT;{Input points and output point  for the interpolation}
  50.  
  51. DEFINE X,Y and Z arrays.
  52. Pick the points for interpolation X0, Y0.
  53. Z0:=INTERP2D(X0,Y0,X,Y,Z,NUMX,MUNY);
  54.  
  55.  
  56.  
  57.  
  58. EXAMPLE
  59.  
  60. The example routine defines Z as Z=sin(x)*cos(y) on a 10 X 10 grid.  X and Y range from 0 to (9/5)ΓêÅ.  The output of ALAG3D is then used to generate a
  61. three dimensional plot of the function from X = 0 to (28/15)ΓêÅ, Y = 0 to (19/20)ΓêÅ.  The function is plotted using the procedure PLOT3.  PLOT3 will scale and plot to a point or scale and move to a point from the current pen position given the point (x,y,z)  and three Euler angles (alpha,beta,gamma). It assumes the screen limits have been defined as XMIN,XMAX,YMIN,YMAX.  To get PLOT3 to move to a point without drawing set color equal to 0.  If color is 1 then the line will be drawn.  The calling procedure for PLOT3 is 
  62. PLOT3(x,y,z,alpha,beta,gama,color);
  63.  
  64. x, y, z, alpha, beta, and gamma are of type float. Color is of type integer.
  65.  
  66.  
  67.  
  68. REFERENCES
  69.  
  70. James, M. L., Smith, G. M., Wolford, J. C.:  Applied Numerical Methods for Digital Computation With Fortran and CSMP, Harper and Row, New York, 1977.
  71.  
  72. Spiegel, M. R..: Mathematical Handbook of Formulas and Tables,  McGraw-Hill, New York, 1968.
  73.  
  74.  
  75. -- part contents for background part 26
  76. ----- text -----
  77. 5
  78.  
  79. -- part contents for background part 6
  80. ----- text -----
  81. PROGRAM TEST_ALAG3D;
  82. TYPE
  83. FLOAT = REAL;
  84. DATA = ARRAY[1..20] OF FLOAT;
  85. RECTANG = ARRAY[1..20, 1..20] OF FLOAT;
  86. FIXXED = INTEGER;
  87. VAR
  88. X0, Y0, Z0, XMAX, XMIN, YMAX, YMIN, A, B, C : FLOAT;
  89. X, Y : DATA;
  90. Z : RECTANG;
  91. NX, NY, K, K2 : FIXXED;
  92. PROCEDURE see;
  93. VAR
  94. R : Rect;
  95. BEGIN
  96. HideAll;
  97. SetRect(R, 0, 35, 550, 400);
  98. SetDRAWINGRect(R);
  99. ShowDRAWING;
  100. END;
  101. FUNCTION alag3d (X0, Y0 : FLOAT;
  102. X, Y : DATA;
  103. Z : RECTANG;
  104. NX, NY : FIXXED) : FLOAT;
  105. VAR
  106. I, J, K, L, M : FIXXED;
  107. U1, U2, U3, U4, V1, V2, V3, V4 : FLOAT;
  108. LA1, LA2, LA3, LA4, LA5, LA6, LA7, LA8 : FLOAT;
  109. QUIT : BOOLEAN;
  110. BEGIN
  111. L := 2;
  112. M := NX - 1;
  113. Quit := FALSE;
  114. REPEAT
  115. REPEAT
  116. I := (L + M) DIV 2;
  117. UNTIL ((L - I) <= 0);
  118. IF (((L - I) = 0) OR ((X[I] - X0) = 0.0)) THEN
  119. Quit := True;
  120. IF (NOT quit) THEN
  121. IF ((X[I] - X0) > 0) THEN
  122. M := I
  123. ELSE
  124. BEGIN
  125. IF (X[I + 1] - X0 < 0) THEN
  126. L := I
  127. ELSE
  128. Quit := True
  129. END;
  130. UNTIL (Quit);
  131. Quit := FALSE;
  132. L := 2;
  133. M := NY - 1;
  134. REPEAT
  135. REPEAT
  136. J := (L + M) DIV 2
  137. UNTIL ((L - J) <= 0);
  138. IF (((L - J) = 0) OR ((Y[J] - Y0) = 0.0)) THEN
  139. Quit := True;
  140. IF (NOT quit) THEN
  141. IF ((Y[J] - Y0) > 0) THEN
  142. M := J
  143. ELSE
  144. BEGIN
  145. IF (Y[J + 1] - Y0 < 0) THEN
  146. L := J
  147. ELSE
  148. Quit := True;
  149. END;
  150. UNTIL (Quit);
  151. U1 := ((X[I - 1] - X[I]) * (X[I - 1] - X[I + 1]) * (X[I - 1] - X[I + 2]));
  152. U1 := (X0 - X[I]) * (X0 - X[I + 1]) * (X0 - X[I + 2]) / U1;
  153. U2 := ((X[I] - X[I - 1]) * (X[I] - X[I + 1]) * (X[I] - X[I + 2]));
  154. U2 := (X0 - X[I - 1]) * (X0 - X[I + 1]) * (X0 - X[I + 2]) / U2;
  155. U3 := ((X[I + 1] - X[I - 1]) * (X[I + 1] - X[I]) * (X[I + 1] - X[I + 2]));
  156. U3 := (X0 - X[I - 1]) * (X0 - X[I]) * (X0 - X[I + 2]) / U3;
  157. U4 := ((X[I + 2] - X[I - 1]) * (X[I + 2] - X[I]) * (X[I + 2] - X[I + 1]));
  158. U4 := (X0 - X[I - 1]) * (X0 - X[I]) * (X0 - X[I + 1]) / U4;
  159. V1 := ((Y[J - 1] - Y[J]) * (Y[J - 1] - Y[J + 1]) * (Y[J - 1] - Y[J + 2]));
  160. V1 := (Y0 - Y[J]) * (Y0 - Y[J + 1]) * (Y0 - Y[J + 2]) / V1;
  161. V2 := ((Y[J] - Y[J - 1]) * (Y[J] - Y[J + 1]) * (Y[J] - Y[J + 2]));
  162. V2 := (Y0 - Y[J - 1]) * (Y0 - Y[J + 1]) * (Y0 - Y[J + 2]) / V2;
  163. V3 := ((Y[J + 1] - Y[J - 1]) * (Y[J + 1] - Y[J]) * (Y[J + 1] - Y[J + 2]));
  164. V3 := (Y0 - Y[J - 1]) * (Y0 - Y[J]) * (Y0 - Y[J + 2]) / V3;
  165. V4 := ((Y[J + 2] - Y[J - 1]) * (Y[J + 2] - Y[J]) * (Y[J + 2] - Y[J + 1]));
  166. V4 := (Y0 - Y[J - 1]) * (Y0 - Y[J]) * (Y0 - Y[J + 1]) / V4;
  167. LA1 := (Z[I - 1, J - 1] * U1 + Z[I, J - 1] * U2 + Z[I + 1, J - 1] * U3 + Z[I + 2, J - 1] * U4) * V1;
  168. LA2 := (Z[I - 1, J] * U1 + Z[I, J] * U2 + Z[I + 1, J] * U3 + Z[I + 2, J] * U4) * V2;
  169. LA3 := (Z[I - 1, J + 1] * U1 + Z[I, J + 1] * U2 + Z[I + 1, J + 1] * U3 + Z[I + 2, J + 1] * U4) * V3;
  170. LA4 := (Z[I - 1, J + 2] * U1 + Z[I, J + 2] * U2 + Z[I + 1, J + 2] * U3 + Z[I + 2, J + 2] * U4) * V4;
  171. LA5 := (Z[I - 1, J - 1] * V1 + Z[I - 1, J] * V2 + Z[I - 1, J + 1] * V3 + Z[I - 1, J + 2] * V4) * U1;
  172. LA6 := (Z[I, J - 1] * V1 + Z[I, J] * V2 + Z[I, J + 1] * V3 + Z[I, J + 2] * V4) * U2;
  173. LA7 := (Z[I + 1, J - 1] * V1 + Z[I + 1, J] * V2 + Z[I + 1, J + 1] * V3 + Z[I + 1, J + 2] * V4) * U3;
  174. LA8 := (Z[I + 2, J - 1] * V1 + Z[I + 2, J] * V2 + Z[I + 2, J + 1] * V3 + Z[I + 2, J + 2] * V4) * U4;
  175. alag3d := 0.5 * (LA1 + LA2 + LA3 + LA4 + LA5 + LA6 + LA7 + LA8);
  176. END;
  177. PROCEDURE plot3 (x, y, z, alpha, beta, gamma : FLOAT;
  178. color : FIXXED);
  179. VAR
  180. x1, x2, y1, y2, x3, y3 : FLOAT;
  181. xint, yint : integer;
  182. BEGIN
  183. x1 := x * cos(gamma) - y * sin(gamma);
  184. y1 := x * sin(gamma) + y * cos(gamma);
  185. y2 := y1 * cos(beta) - z * sin(beta);
  186. x3 := x1 * cos(alpha) - y2 * sin(alpha);
  187. y3 := x1 * sin(alpha) + y2 * cos(alpha);
  188. yint := 200 - round((y3 - ymin) * 200 / (ymax - ymin));
  189. xint := round((x3 - xmin) * 495 / (xmax - xmin));
  190. IF (color = 1) THEN
  191. lineto(xint, yint)
  192. ELSE
  193. moveto(xint, yint);
  194. END;
  195. BEGIN
  196. see;
  197. XMAX := 6.0;
  198. XMIN := -6.0;
  199. YMAX := 7.0;
  200. YMIN := -7.0;
  201. FOR K := 1 TO 10 DO
  202. BEGIN
  203. X[K] := 3.1415926 * (K - 1.0) / 5.0;
  204. FOR K2 := 1 TO 10 DO
  205. BEGIN
  206. Y[K2] := 3.1415926 * (K2 - 1.0) / 5.0;
  207. Z[K, K2] := SIN(X[K]) * COS(Y[K2]);
  208. END;
  209. END;
  210. A := 3.14;
  211. B := 0.707;
  212. C := 0.707;
  213. FOR K := 1 TO 15 DO
  214. BEGIN
  215. X0 := 3.1415926 * (K - 1.0) / 7.5;
  216. FOR K2 := 1 TO 20 DO
  217. BEGIN
  218. Y0 := 3.1415926 * (K2 - 1.0) / 10;
  219. Z0 := ALAG3D(X0, Y0, X, Y, Z, 10, 10);
  220. IF (K2 = 1) THEN
  221. PLOT3(X0, Y0, Z0, A, B, C, 0)
  222. ELSE
  223. PLOT3(X0, Y0, Z0, A, B, C, 1);
  224. END;
  225. END;
  226. FOR K := 1 TO 3 DO
  227. BEGIN
  228. Y0 := 3.1415926 * (K - 1.0) * 0.95;
  229. FOR K2 := 1 TO 40 DO
  230. BEGIN
  231. X0 := 3.1415926 * (K2 - 1.0) * (28.0 / 15.0) / 39;
  232. Z0 := ALAG3D(X0, Y0, X, Y, Z, 10, 10);
  233. IF (K2 = 1) THEN
  234. PLOT3(X0, Y0, Z0, A, B, C, 0)
  235. ELSE
  236. PLOT3(X0, Y0, Z0, A, B, C, 1);
  237. END;
  238. END;
  239. END.